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Abstract 

If the dynamics of an evolutionary differential equation system possess a low- 
dimensional, attracting, slow manifold, there are many advantages to using this 
manifold to perform computations for long term dynamics, locating features such as 
stationary points, limit cycles, or bifurcations. Approximating the slow manifold, 
however, may be computationally as challenging as the original problem. If the 
system is defined by a legacy simulation code or a microscopic simulator, it may 
be impossible to perform the manipulations needed to directly approximate the 
slow manifold. In this paper we demonstrate that with the knowledge only of a 
set of "slow" variables that can be used to parameterize the slow manifold, we can 
conveniently compute, using a legacy simulator, on a nearby manifold. Forward 
and reverse integration, as well as the location of fixed points are illustrated for a 
discretization of the Chafee- Infante PDE for parameter values for which an Inertial 
Manifold is known to exist, and can be used to validate the computational results. 
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1 Introduction 

Certain dissipative evolutionary equations possess low-dimensional, attracting invariant 
manifolds which govern their long term dynamics. Such a manifold is readily apparent 
for a system given in the singularly perturbed form: 



dy/dt = y' = f(y,z), f G 5T +m h-> 3? 



dz/dt = z' = g(y, z)/e, g G 5T +m i-> K m (1) 

where / and g are such that, for small positive e, z is rapidly attracted to the region 
z' = 0(1) and dg/dz is non-singular. Since z = z {y) + 0(e) where z (y) is the solution 
of g(y, z (y)) = 0, the slow manifold is given by the solution of 

y' = f(y,z (y)) + O(e) (2) 

so that we can "easily" compute an approximation to it. The one-dimensional slow 
manifold is parameterized here by y (other parameterizations are also possible). 

In more complicated cases, approximations to the slow manifold may not be so appar- 
ent; yet within such manifolds the system dynamics can still be described by a lower-order 
differential equation - the reduced system. Methods for approximating such manifolds 
have been the subject of intense research in communities ranging from reactive flow mod- 
eling (e.g. [SIIS1IZI) to inertial manifolds for dissipative PDEs (e.g. [HI El II])- If we are 
able to somehow constrain the dynamics to a slow manifold, stable numerical integration 
could be performed with larger stepsizes than would be possible in the original system. 
Furthermore, many global properties of the original system are (approximately) inherited 
by the reduced system; these include stationary points, limit cycles, and bifurcations and 
may be computable more easily on the slow manifold. Unfortunately, approximating the 
slow manifold may be as computationally challenging as the original problem. 

In our work we seek an approximation to such a manifold that is (a) simple to obtain 
on the fly during numerical computations, and (b) only requires evaluations of time 
derivatives of the state, such as would be available from a legacy code. Our starting 
point is the assumption that, given a basis for the full set of variables in the problem, a 
subset of this basis can be used to parameterize the slow manifold and our approximation 
of it, as y did in our example above. In some applications, such as when the full system is 
described by a microscopic simulation, the subset used to parameterize the slow manifold 
might be called "macroscopic observables" ; such observables could be the pressure field 
in kinetic theory based flow simulation, or a concentration field in the kinetic Monte 
Carlo simulation of a chemical reaction. 

We may start with a finite-dimensional system (an ODE) or an infinite-dimensional 
system (a PDE, for example). In the latter case we will have to introduce a finite- 
dimensional approximation before commencing computation (in effect, the "method of 
lines"). Suppose that the system is represented by (approximated by) the system 

u t = L(u) (3) 

where u = {ui}, i = 1, ■ • • , N. For ease of presentation, let us assume that the equation 
has already been transformed to a suitable basis so that U = {Ui — i — 1, • • • , n} 
parameterizes the slow manifold. In some sense we are assuming that the remaining 
variables, V = {Vi = u n+ i, i = 1, ■ • ■ , iV — n = m}, are the "fast" ones that are quickly 
"slaved" to U; we will return to this assumption. The split system is 

1e =Li{u > v) (4) 

dV 

^ = L 2 (U,V). (5) 

We now view the solutions of the system as families of curves in the U — V state space, 
as illustrated in Figure [T] - although such figures are potentially misleading because we 
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have to remember that each axis represents a multi-dimensional space. The essential 
feature of the figure is that one member of the family of solutions is a "slow" manifold 
with no high- curvature region, while other members of the family of solutions approach 
this slow manifold relatively rapidly. 

Since the slow manifold, M, can be parameterized by the slow variables, U, points 
on M, (U, Vm(U)), must be uniquely determined by U - that is, the curve cannot "fold" 
in the region of interest. If we had a scheme for approximating the value of V — Vm(U) 
for each U (as we did in the singularly perturbed example above) we could, for example, 
apply a numerical integration method to just the U variables, computing the equivalent 
values of Vm{U) only as needed by the integration scheme. This is the main point of 
our approach: one does not compute the entire manifold a priori, but only computes 
it pointwise, "on demand" as required by the low- dimensional integration code (or by 
algorithms performing other numerical tasks, such as fixed point computation). 

From the assumptions, we suspect that an approximation to the slow manifold can 
be found by computing the value of V — Va(U) (as shown in Figure EJ) for which the 
time derivatives of the V components are zero. Here we will compute on this "steady 
manifold"; for the appropriate basis choice this steady AIM is not too far from the 
slow manifold. Better approximations, based on higher order expansions of singularly 
perturbed equations, can also be used in a legacy code context, and will be the subject 
of further work. 

We are particularly interested in the case where we do not have "access" to the 
differential equations directly because, for example, we have a legacy simulation code, 
or the system is the unknown closure of a microscopic simulation model (kinetic Monte 
Carlo, molecular dynamics). Then, the only computational possibility we have is to 
integrate the full system for a short time in what we call an inner integration step. In 
this case, we can define an AIM by requiring that the chord of the inner step has zero 
change in the V components, that is, it is "horizontal." This can be computed iteratively 
by performing an inner integration over a small step of length h and then projecting the 
solution back to the specified value of U (by simply resetting the values of U to their 
values at the start of the step) and repeating. If h is small enough, this iteration will 
converge if the solution family in the neighborhood of (U, Va{U)) is locally attracted to 
the solution that passes through (U, Va(U)). This can be seen by noting that one step 
of the process performs the mapping 

V m+1 = V m + gL 2 (U, V m ) + 0(h 2 ) (6) 

so we have 

V m+1 -V m = {I + h-^)(V m - Kn-i) + 0(h 2 ) 

which implies convergence for small enough h if the eigenvalues of dLi/dV are in the 
negative half plane. This property can form the basis of alternative algorithms to approx- 
imate the "steady AIM:" matrix free fixed point algorithms, like the Recursive Projection 
Method or GMRES ^01 HTj can be applied to accelerate the computation of the fixed 
point of eq. (JHJ). 

In the following sections we will illustrate the use of this technique on the Chafee- 
Infante reaction-diffusion equation 

du d 2 u n 
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Figure 1: Family of Solutions and Slow Manifold Schematic 
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Figure 2: Approximate Inertial Manifold Definition 
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with v = (2.5)~ 2 and u(Q,t) = u(tt, t) = 0. This is known to possess an inertial manifold 
of dimension two (in effect, the two-dimensional unstable manifold of the origin, and its 
closure). Although we know the differential equations in this example, we are not going 
to make explicit use of the knowledge in our computational method. We will only use it 
as if we were given a legacy code for evaluating time derivatives. 

We first discretize the equations in space. Since we are not interested in the issue 
of the best spatial discretization, we use simple finite difference methods over N equally 
spaced points, so that the variables are Ui(t) = u(xi,t),i = 1, • • • , N, where Xi = m/(N + 
1). These variables are chosen for convenience in the calculation. The resulting ODEs 
are the usual: 

-T[ = -^j^{ui-i - 2ui + u i+1 ) + Ui - i = l,---,N (8) 

where uo = un+i = 0. If we had a legacy code or a microscopic model, the variables 
would be the ones that happened to come with the code or model. In this example, no 
subset of the {u{\ variables is suitable for defining the AIM (since the slow manifold 
varies rapidly as a function of each Uj) so we will use an "observation basis" in which a 
linear combination of the variables will parameterize the AIM. In this case, we can use 
a basis formed by sin(mx), m = 1, • • • , N. (These are the unnormalized eigenvectors of 
eq. (JBJ) when u = 0.) The modified variables are 



N 



~>ijUj 



(9) 



where fa = is the basis given by faj = sin(ixj). The first two a\ can parameterize 

the slow manifold, and it is not necessary to calculate the rest. 

We now present a technique for approximating d(ai, GS2) / 'dt on the slow manifold 
given (ai,a 2 ). This approximation can then be used to implement time integration, 
stability analysis, or other numerical procedures on the system constrained to the AIM. 
The general method consists of 

1. Start with a prescribed value of (ai, g^). 

2. Compute the values of such that eq. (JHJ) is satisfied and the local derivative 
of the full {{ui} system) is "horizontal" in the other components of the basis (in 
this example, 03, • • • , ajv)- This can be done in a number of ways: 

(a) Use eq. (JSJ) and eq. to compute dai/dt, 3 < i =< N and then solve for the 
values of aj that makes these zero using Newton iteration. This can be done 
directly when the equations are available (it is done in the example illustrated 
here), or can be done through matrix- free based contraction mappings if the 
equations are not explicitly available. 

(b) If we only have a legacy code or a microscopic simulator of the full system, 
use iteration eq. © repeatedly to find the values of aj, 3 < j < N such that 
the chords of those aj are zero. 

(c) Conceptually (since this is not practical for a legacy code) one could imple- 
ment a Lagrange multiplier, evolving the dynamics while constraining the 
projection of the solution on (ai,a 2 ). This is reminiscent of techniques like 
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SHAKE used to "prepare" molecular dynamics simulations The ap- 

proach described immediately above is a way of effectively implementing what 
amounts to such a Lagrange multiplier constrained integration to a legacy 
simulator. 

3. Compute the derivatives (or chord slope) of the full system from the given 
values of a\, a 2 and the now computed a,, 3 < j < N (actually they have probably 
been computed in the previous step). 

4. Compute the "(01,02)" components of the derivative by applying eq. to the 
{ui} derivatives. These are the approximations to the time derivatives of (01,02) 
on the steady AIM. 

In the next section we will use this technique to integrate eq. (jHJ) both forward 
and backward in time on our two-dimensional steady AIM, and compare it with the 
integration of the full system and in the subsequent section we will use it to compute the 
steady states directly by performing a Newton iteration on the two-dimensional steady 
AIM. 

2 Integration on an AIM of the Reaction-Diffusion 
Equation 

In [T] we introduced projective integration which uses computation of the chord slope 
obtained by integration of a legacy code or of a microscopic model in place of derivatives 
for performing large projective integration steps on the slow components. If we were 
working with legacy codes or microscopic simulators, we would use that technique in our 
"on manifold" integrations. However, we have chosen an example for which we know the 
equations of the detailed system, so that we can compare the "true" integration of the 
system with the approximation on the steady AIM we have defined. 

We integrated the full system eq. (JHJ) using the automatic Runge-Kutta method with 
the Dormand-Prince pair of formulae known as RK45 and available in MATLAB as 
ode45. We also used the same method to integrate just (ai,o 2 ) on the AIM, using the 
technique described in the previous section to approximate the derivatives of (ai, 02). It 
is possibly better to view this as the (approximate) integration of the projection on the 
(01,02) "observables plane" of the true dynamics. The results for (01,02) are shown in 
Figure El The integration was started from six different points near the origin in the 
(ai, a 2 ) plane (which is an unstable steady state). All but one approach the stable steady 
state at approximately (1.12,0) but the one that starts at (0,0.05) stays on the invariant 
submanifold ai = and moves to the saddle point at about (0,0.7). Since the origin is 
on the inertial manifold of eq. (jHJ), the starting points are also very close to it so that 
the RK45 integration of the full system gives a good approximation of the solutions on 
the inertial manifold and provides a picture of the manifold itself. In Figure El the solid 
line is the RK45 solution of the full system, while the dashed line is the integration on 
our steady AIM. 

The full system in eq. (JB1) is rapidly damped in its fastest components, and so it 
would not be feasible to numerically integrate it in the reverse time direction. However, 
the (01,02) differential system on the AIM does not have these fast components, so 
it can be integrated "backwards." The dotted arcs in Figure El are the results of a 
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reverse integration in the (ai, <22)-plane starting from a point on the RK45 solution of 
the full system shortly before the stationary point is reached (one can't start too close 
to the stationary point because the trajectory chosen would be too sensitive to small 
perturbations). 

As we can see in FigureEl the forward and reverse solutions on the AIM are fairly good 
approximations to the components of the "true" solution on the IM. We do not expect 
the other components, aj, i > 3 to be good approximations. This is shown in Figure 0] 
which shows the values of 03 plotted against a\ for each of the trajectories. (The reverse 
integration trajectories are almost indistinguishable from the forward trajectories for the 
(ai,a 2 ) integration.) 



3 Steady State Computation on the AIM 

The procedure described in Section 1 computes da/dt = p(a) where a = (ai,a 2 ) is the 
parameterization of the slow manifold. Using any standard techniques we can look for 
zeros of p(a) to identify stationary states of the system. Since the dimension is low, 
we can use Newton's method, computing the approximate partial derivatives by finite 
differencing. Many better methods exist, but our purpose here is simply to show that 
the reduced system can be used directly in any conventional numerical process. 

Tabled shows the sequence of iterates for Newton's method starting at three different 
point in the (ai, a 2 )-plane and iterating until changes were less than 10 -5 in the LI norm. 
The eigenvalues of J = dp/ da as the iteration proceeds are also shown. For comparison, 
the two leading eigenvalues of the full system eq. (JBJ) at steady state for each of the three 
cases are (0.8403,0.3648), (0.2204,-0.7118), and (-1.4491,-1.5392) respectively. As can 
be seen, the three stationary points, (0,0), (0,0.7056), and (1.1206,0) are source, saddle, 
and stable (sink), respectively. Note that the stationary states on the steady AIM are 
necessarily on the slow manifold since all of the derivatives are zero at these points. 



4 Discussion 

We have demonstrated that it is possible to perform low-dimensional (macroscopic) com- 
putations on an AIM (more precisely, on observations of an AIM) based on choosing a 
suitable parameterization of the low dimensional slow manifold. That parameterization 
must be chosen so that the slow manifold does not "fold over it" . It must also be chosen 
so that the induced AIM is reasonably close to the true slow manifold. In the example 
we discussed, the first two eigenvectors of the linearization of the problem at a particular 
solution value (the origin) were chosen to parameterize the manifold, since at that solu- 
tion value they are tangent to the true slow manifold. As long as the solution does not 
stray too far from that region (compared to the non-linearities present) these directions 
provide a reasonable parameterization to the slow manifold elsewhere. The steady AIM 
for these two variables is illustrated in Figure El which plots 03 against a\ and a 2 on the 
AIM. This AIM is a reasonable approximation of the slow manifold, (which is shown is 
Figure 1 of |13j). It is clear from this figure that choosing, say, a 2 and 03 to characterize 
the slow manifold would have been bad since ai(a 2 , a 3 ) is multivalued for some a 3 . (See 
the line a 2 = l,a 3 = —0.1 that is indicated in the figure. It intersects the AIM twice 
within the region plotted.) This is true even near the origin, where this AIM is a good 
approximation to the inertial manifold. 
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Table 1: Sequence of Newton Iterates 



Case 1 


a\ 


a 2 


Ai 


A 2 


0.2000 


0.2000 


0.7171 


0.1971 


-0.0838 


-0.1946 


0.7737 


0.2684 


0.0159 


0.0551 


0.8351 


0.3571 


-0.0002 


-0.0009 


0.8403 


0.3648 


0.0000 


0.0000 


0.8403 


0.3648 


Case 2 


a-y 


a 2 


Ai 


A 2 


0.1000 


0.7500 


0.1666 


-0.8904 


-0.0405 


0.7271 


0.2069 


-0.7902 


0.0063 


0.7087 


0.2359 


-0.7293 


-0.0001 


0.7057 


0.2405 


-0.7203 


0.0000 


0.7056 


0.2407 


-0.7200 


Case 3 


a\ 


a 2 


Ai 


A 2 


1.0000 


0.1000 


-0.8687 


-1.5201 


1.1707 


-0.0462 


-2.0297 


-1.6832 


1.1255 


-0.0042 


-1.6847 


-1.6599 


1.1206 


-0.0000 


-1.6529 


-1.6529 


1.1206 


-0.0000 


-1.6528 


-1.6526 
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The focus of this work was on the use of a legacy simulator to approximate the slow 
manifold on the fly, as dictated by the needs of numerical analysis tools employed for 
computations on it. The local nature of the approximation should be contrasted to "off 
line" algorithms that attempt to approximate the entire manifold first (see the extensive 
discussion in ^3] as well as J2|). Here we pursued the simplest approximate manifold 
one can find by constraining a legacy code. An important issue that was only tangentially 
mentioned here was the selection of good basis functions (or macroscopic observables in 
the case of atomistic inner simulators) that parameterize the manifold; statistical data 
analysis techniques have an important role to play in this. Better algorithms, resulting 
from the implementation of higher order approximations to the slow manifold (requiring 
a higher order derivative to vanish) in a legacy code context are currently being explored. 

Acknowledgements This work was partially supported by an NSF/ITR grant and 
by the AFOSR Dynamics and Control (Dr. B. King). 
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Figure 3: Trajectories in (ai, a 2 )-plane. Solid: Inertial manifold, dashed - AIM forward, 
dotted - AIM reverse. 
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